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We present the first two-loop calculation of the heavy quark energy shift in lattice nonrelativistic 
QCD (NRQCD). This calculation allow us to extract a preliminary prediction of m^irni,, rif = 
5) = 4.25(12) GeV for the mass of the b quark from lattice NRQCD simulations performed with 
a lattice of spacing a = 0.12fm. Our result is an improvement on a previous determination of 
the b quark mass from unquenched lattice NRQCD simulations, which was limited by the use 
of one-loop expressions for the energy shift. Our value is in good agreement with recent results 
of m.b(mb) — 4.163(16) GeV from QCD sum rules and m&(m&, »/ = 5) = 4.165(23) GeV from 
realistic lattice simulations using highly-improved staggered quarks. We employ a mixed strategy 
to simplify our calculation. Ghost, gluon and counterterm contributions to the energy shift and 
mass renormalisation are extracted from quenched high-beta simulations whilst fermionic contri- 
butions are calculated using automated lattice perturbation theory. Our results demonstrate the 
effectiveness of such a strategy. 
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1. Introduction 

The precise theoretical and experimental determination of quark masses is an important com- 
ponent of high-precision tests of the Standard Model of particle physics. One current focus for tests 
of the Standard Model is the unitarity of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, which 
describes flavour-changing quark transitions. Quark masses serve as an input into the tests of CKM 
matrix unitarity; the mass of the b quark is used in the extraction of the CKM matrix element \V u b\ 
from inclusive semileptonic decays of B mesons [JTJ] _ 

Recent high-precision calculations of the b quark mass using realistic lattice QCD simulations 
[0J and perturbative QCD combined with experimental results [||] are in good agreement, obtaining 
values of nib(mb, rif = 5) = 4.165(23) GeV and mb(mb) = 4.163(16) GeV respectively. The lat- 
tice result was obtained using the highly improved staggered quark (HISQ) action for the valence 
quarks, with the improved relativistic Asqtad action used for the sea quarks. HISQ is a highly cor- 
rected version of the standard staggered quark action that retains a chiral symmetry on the lattice 
[Q]. Most current lattice studies of b quarks use an effective field theory, such as nonrelativistic 
QCD (NRQCD), for the valence heavy quarks. Simulating both valence and sea quarks with rela- 
tivistic actions allows a much greater precision, but is only now becoming possible with the advent 
of finer lattices and highly improved actions. However, even on the very finest lattices with HISQ 
heavy quarks an extrapolation to the heavy quark mass is still required [0]. 

Our calculation improves on a previous determination of m/,(m^) = 4.4(3) GeV from un- 
quenched lattice QCD simulations using NRQCD valence b quarks The dominant error in 
that calculation arose from the use of one-loop perturbation theory in the matching between lat- 
tice quantities and the continuum result. By introducing a mixed strategy incorporating high-beta 
quenched simulations and automated lattice perturbation theory, we perform the first ever such two- 
loop calculation in NRQCD. This serves a two-fold purpose. Firstly our calculation demonstrates 
the effectiveness of employing an efficient mixed perturbation theory/high-beta simulation method 
for higher order perturbative quantities. Secondly our result allows us to obtain a more precise 
prediction for the b quark mass from lattice NRQCD simulations. 

1.1 Heavy quarks on the lattice 

Currently available lattices are too coarse to directly simulate b quarks, because the Compton 
wavelength of the b quark is smaller than the lattice spacing. One common approach to solving 
this problem is to introduce a nonrelativistic effective action, NRQCD, for which the discretization 
errors are under control and which can be systematically improved by including extra operators. 

NRQCD is constructed by integrating out dynamics at the scale of the heavy quark mass and 
then using the Foldy-Wouthuysen-Tani transformation to write the action as an expansion in the 
inverse heavy quark mass [|6||. We use an NRQCD action correct to <^(l/m 2 , v 4 ), where v is the 
relative internal velocity of the bound-state heavy quarks. A detailed derivation of the action we 
use is given in |{7|]. The lattice NRQCD action can be written 

Snrqcd = £ Y + 0, *) [V( X > T ) ~ Vv*, T ~ ^ > ( U ) 
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with 




(1.2) 



Here the leading nonrelativistic kinetic energy is Hq = —A^/2M. The correction term 8H con- 
tains higher order terms in the 1 /M expansion: the improved chromoelectric and chromomagnetic 
interactions, the leading relativistic kinetic energy correction and discretization error corrections. 
The integer n is introduced as a stability parameter. 



2. Calculating the b quark mass 

Quark confinement ensures that quark masses are not physically measurable quantities, so the 
notion of quark mass is a theoretical construction. A wide range of quark mass definitions ex- 
ist, often tailored to exploit the physics of each particular process. One common choice of quark 
mass is the pole mass, defined as the pole in the renormalized heavy quark propagator. However, 
the pole mass is a purely perturbative concept and suffers from infrared renormalon ambiguities 
To avoid these ambiguities, experimental results are usually quoted in the modified Min- 
imal Subtraction (MS) scheme, which is renormalon ambiguity free. Lattice calculations use the 
renormalon-free bare lattice mass. These different quark mass definitions must be matched to en- 
able meaningful comparison. We match bare lattice quantities to those in the MS scheme using 
the pole mass as an intermediate step. Any renormalon ambiguities cancel in the full matching 
procedure between the lattice quantities and the MS mass. We extract the MS mass from lattice 
simulation data in a two-stage process. We first relate lattice quantities to the pole mass and then 
match the pole mass to the MS mass evaluated at a scale equal to the b quark mass. 

2.1 Extracting the pole mass 

We determine the pole mass using two independent methods. The first method relates the pole 
mass, Mf ole , to the experimental T mass, 

M ex P t = 9 46030(26) GeV [jUjl, using the heavy quark 

energy shift, Eq\ 

2M p b ole = M™ pt - (£ sim (0) - 2E ). (2. 1) 

Here £" sim (0) is the energy of the T meson at zero momentum, extracted from lattice NRQCD 
simulations. The quantity (£ sim (0) — 2Eq) corresponds to the 'binding energy' of the meson in 
NRQCD. We use a value of £ sim = 0.515(3) GeV, obtained from a lattice NRQCD simulation run 
by the HPQCD collaboration on a "coarse" MILC ensemble, with lattice spacing a = 1.647(3) 



GeV 1 [jl2j. For further details of the configuration ensemble see [|13(, |14J 



The second method directly matches the pole mass to the bare lattice mass in physical units, 

M latt 

(a), via the heavy quark mass renormalisation, Zj^ tt , 

M pole = Z^\pia,M^(a))M^{a). (2.2) 

We employ a mixed strategy to calculate Eq and Zj^ tt perturbatively. The fermionic contribu- 
tions to Eq, shown on the left-hand side of Figure [IJ are calculated using two-loop automated lattice 
perturbation theory. All other contributions, shown on the right-hand side of Figure [j], are extracted 
from high-beta quenched simulations. 
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Figure 1: Contributions to Eq and Z]|". The four fermionic contributions calculated using automated lattice 
perturbation theory are shown on the left. The diagrams on the right are extracted from high-j3 simulation. 
Blue lines are heavy quarks, green are gluons and red are sea quarks. Large brown blobs represent the 5 
gluon self energy diagrams and crosses are counterterms. Feynman diagrams reproduced from jl5|]. 



Results were obtained using the NRQCD action of Equation [Tj] for the heavy valence quark, 
HISQ light quarks and the Liischer-Weisz action for the gluons flTrjj , Oj. We used a heavy quark 
mass in lattice units of Ma = 2.8, with a stability parameter of n = 2. 

2.1.1 Automated lattice perturbation theory 

Feynman rules for the NRQCD and HISQ actions are too complicated to be viably derived 
by hand and the resulting Feynman integrals can only be evaluated numerically. We therefore use 
automated lattice perturbation theory, employing HiPPy to derive the Feynman rules and HPsrc 



to evaluate the four diagrams [18, 19]. To control the highly-peaked IR behaviour of the Feynman 
integrands, we introduce a gluon mass. Although in general a non-zero gluon mass cannot be 
used in calculations beyond one-loop, this issue concerns only diagrams containing ghost-gluon 
vertices. In our calculation, these diagrams are handled by the high-beta simulation, allowing us 
to use a gluon mass for the fermionic contributions. The light quark diagrams in Figure [j] were 
calculated using five different light quark masses and extrapolated to zero light quark mass. We 
verified that the appropriate Ward identity for the 1-loop gluon self-energy was satisfied. 

2.1.2 High-beta simulations 

We perform quenched simulations on l? x T lattices with temporal extent T = 3L, for L = 3 
to L = 10 and twisted boundary conditions to reduce finite size effects and tunnelling between 



QCD vacua [17]. We generate ensembles of configurations for 17 values of j8 from j8 = 9 to j8 = 
120. Since the Green function is not gauge-invariant, we fix to Coulomb gauge using a conjugate 
gradient method. To extract the energy shift and mass renormalisation, we use an exponential fit to 
the heavy quark Green function parametrized as 



G(p,t) =Z v ,exp 



E0 + 2^ + - 



(2.3) 
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where the ellipsis stands for higher order terms that are included in the fits. 

All operators in the NRQCD action are expressed in terms of gauge-covariant Wilson paths 
generated using PYTHON, which greatly enables flexibility and reduces programming errors. The 
heavy quark source is classified in the flavour-smell basis appropriate to the twisted boundary 
conditions. We implement the boundary conditions using a gauge-twist mask whenever a path 
in an operator crosses any spatial boundary. By applying an extra U(l) phase in the mask, we 
can assign an arbitrarily small momentum to the source, enabling both Eq and Z M tt to be reliably 
extracted as a function (j3,L). 

We convert j3 to ay and perform a joint fit to extract the 1- and 2-loop coefficients in the 
L — > oo limit. Simulations were run including tadpole improvement, which significantly reduces 
the magnitude of both 1- and 2-loop coefficients. Results for Z l M tt are good but this work is still in 
progress and we concentrate here on those for Eq. For even L, in Table [I] we compare the tadpole- 
improved 1-loop coefficient from an unconstrained fit to the simulation data for Eq with the exact 
result from automated perturbation theory. To extract the 2-loop coefficient we constrained the 1- 
loop coefficient to be the exact value, but Table |l| shows that the simulation reliably reproduces the 
1-loop results. The number of independent configurations for each (j8,L) was about 300, which we 
can easily increase by 10-fold or more, allowing for much more accurate results at the next stage. 



L 


4 


6 


8 


10 


OO 


psim. 


0.5295(16) 


0.5988(16) 


0.6369(12) 


0.6560(11) 


0.7380(63) 


pth. 


0.5312 


0.6020 


0.6362 


0.6565 


0.7348(3) 



Table 1: Comparison of an unconstrained fit from simulation for the perturbative 1-loop coefficient with 
the automated perturbative calculation. There is no error on the theory calculation as it was done by mode 
summation. The error on the theory extrapolation to L — °° is estimated from a fit. 



2.2 Matching the pole mass to the MS mass 

Although the pole mass is plagued by renormalon ambiguities, these ambiguities cancel when 
lattice quantities are related to the MS mass. This renormalon cancellation is evident in the direct 
matching of the bare lattice mass to the MS mass, 

M m (n) = Z^"( At a,M]f tt (a))Z co i t ( At ,M pole )Mr (a), (2.4) 

as both M MS and M l b att are renormalon-free. The continuum matching parameter, Z^ nt , relates the 
pole mass to the MS mass and has been determined to ff(cc^) [|20||. 

To see that renormalon ambiguities also cancel in when determining the pole mass from the 



energy shift, we equate Equations 2.1 and 2.2 and rearrange them to obtain 



2(z]f M^>) " E ) = M^ xpt - £ sim (0). (2.5) 

The two quantities on the right hand side of the equation are renormalon ambiguity free: M^ pt 
is a physical quantity and £' sim (0) is determined nonperturbatively from lattice simulations. Any 
renormalon ambiguities in the two power series, Zj^ tt and Eq, on the left-hand side of the equation 
must therefore cancel. 
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3. Results 

For the fermionic and quenched contributions to the two-loop heavy quark energy shift we find 
E = 0.7348(3)a v (9*/a) + (1.37(6) -0.0023(l)n f )c%(q*/a) + 0(<4) ■ (3-1) 
We express our result in the V-scheme at a scale q* /a = 3.33, a value determined using the BLM 



procedure in Q21|]. The uncertainties quoted for the one-loop coefficient and the quenched contribu- 
tion to the two-loop coefficient arise from the multi-polynomial fit. For the fermionic contribution 
to the two-loop coefficient, the quoted uncertainty is the statistical error in the numerical evalua- 
tion of the Feynman diagrams. We estimated the coefficient of the & (a s f) term from the quenched 
simulation fits as ~ 1.0(5). 



Inserting this result for the heavy quark energy shift into Equation [2JJ leads to our first pre- 
liminary determination of the b quark mass: 

m ms ( m ms\ = 4 25(12) GeV. (3.2) 



The error is an estimate of ^(ce 3 ) contributions, which dominate the uncertainty in our result. 
Uncertainties arising from systematic and statistical errors in the lattice results, i? sim (0) and Eq, are 
< 1%. We are unable to estimate the systematic error due to 0{a 2 ) artifacts as we have not yet 
finished the calculation for smaller values of a; this work is in progress and entails working with 
different values of Ma in NRQCD. It should be noted that we used a value of Z? sim (0) that was 



generated from lattice NRQCD simulations using the action of Equation |1.1| , but with n = 4. From 
1-loop calculations we estimate the errors associated with this mismatch to be much smaller than 
the dominant &((X^) error. However, this discrepancy will be corrected in future work. 



4. Conclusion 

We have calculated the two-loop heavy quark energy shift in highly-improved NRQCD using 
a mixed approach combining quenched high-beta simulations with lattice perturbation theory. This 
is the first determination of any heavy quark parameter beyond first-order perturbation theory in 
NRQCD, and demonstrates that we are able to extract a more precise prediction of the b quark mass 
from lattice NRQCD simulations than has been previously achieved. Work is currently underway to 
complete our calculation of the mass renormalisation, Zj^ tt , and to extend our results to incorporate 
different heavy quark masses to enable us to effectively estimate iff (a 2 ) uncertainties. We also plan 
to increase the size of the ensembles used in the high-beta analysis by a significant factor. We 
expect these developments will improve further the precision of our result for the b quark mass. 



Acknowledgements 

We would like to thank Christine Davies and Iain Kendall for providing HPQCD simulation 

data. 

We thank the DEISA Consortium (www.deisa.eu), funded through the EUFP7 project RI- 
222919, for support within the DEISA Extreme Computing Initiative. This work has made use of 
the resources provided by the Cambridge High Performance Computing service supported in part 
by the Science and Technology Facilities Council under grant ST/H008861/1. 



6 



The b quark mass from lattice nonrelativistic QCD 



Christopher J. Monahan 



References 

[1] E. Barberio, Inclusive semileptonic B decays, (2006) [arXiv : hep-ex/ 6050 98] 

[2] C. McNeile et al. (HPQCD), High-precision c and b masses, and QCD coupling from current- current 
correlators in lattice and continuum QCD, Phys. Rev. D 82 (2010) 034512 [arXiv: 1004.4285] 

[3] K. G. Chetyrkin et al, Charm and bottom quark masses: an update, Phys. Rev. D 80 (2009) 074010 

[arXiv: 0907 . 2110 [hep-ph]] 

[4] E. Follana et al. (HPQCD), Highly improved staggered quarks on the lattice with applications to 
charm physics, Phys. Rev. D 75 (2007) 054502 [arXiv : hep-lat / 0610092] 

[5] A. Gray et al. (HPQCD & UKQCD), The T spectrum and m h from full lattice QCD, Phys. Rev. D 72 
(2005) 094507 [arXiv : hep-lat / 0507013] 

[6] G. P. Lepage et al., Improved nonrelativistic QCD for heavy-quark physics, Phys. Rev. D 46 (1992) 
4052 

[7] R. R. Horgan et al., Moving nonrelativistic QCD for heavy-to-light form factors on the lattice, 
Phys. Rev. D 80 (2009) 074505 [arXiv: 0906 . 094 5] 

[8] G. T. Bodwin, Y-Q. Chen, Renormalon ambiguities in NRQCD operator matrix elements, 
Phys. Rev. D 60 (1999) 054008 [arXiv : hep-ph/ 9807 4 92] 

[9] 1. 1. Bigi et al., Pole mass of the heavy quark: perturbation theory and beyond, Phys. Rev. D 50 
(1994) 2234 [arXiv : hep-ph/ 94 02 3 60] 

[10] M. Beneke, V.M. Braun, Heavy quark effective theory beyond perturbation theory: renormalons, the 
pole mass and the residual mass term, Nucl. Phys. B 426 (1994) 301 [arXiv : hep-ph/ 94 02 3 64] 

[11] K. Nakamura et al. (Particle Data Group), Review of Particle Physics, J. Phys. G 37 (2010) 075021 

[12] C. T. H. Davies, Private communication, 28/09/2010 

[13] A. Bazavov et al., Full nonperturbative QCD simulations with 2+1 flavors of improved staggered 
quarks, Rev. Mod. Phys. 82 (2010) 1249 [arXiv: 903 .35 98] 

[14] C. T. H. Davies et al. (HPQCD), Precise determination of the lattice spacing in full lattice QCD, 
Phys. Rev. D 81 (2010) 034506 [arXiv: 910 . 12 2 9] 

[15] Q. Mason, H. D. Trottier, R. R. Horgan, High precision fundamental constants using lattice 

perturbation theory, in proceedings of XXIIrd International Symposium on Lattice Field Theory, PoS 

(LAT2005) Oil 

[16] Zh. Hao et al, Unquenching effects on the coefficients of the Liischer-Weisz action, Phys. Rev. D 76 
(2007) 034507 [arXiv : 7 5 . 4 6 6 0] 

[17] M. Luscher, P.Weisz, Efficient numerical techniques for perturbative lattice gauge theory 
calculations, Nucl. Phys. B 266 (1986) 309 

[18] A. Hart et al, Automated generation of lattice QCD Feynman rules Comp. Phys. Comm. 180 (2009) 
2698 [arXiv: 0904 . 0375] 

[19] T. C. Hammant et al, Improved automated lattice perturbation theory in background field gauge, in 
proceedings of XXVIII International Symposium on Lattice Field Theory 

[20] K. Melnikov, T. van Ritbergen, The three-loop relation between the MS-bar and the pole quark 
masses, Phys. Lett. B 482 (2000) 99 [arXiv : hep-ph/ 9 912 3 91] 

[21] K. Y. Wong, H. D. Trottier, R. M. Woloshyn, Perturbative Wilson loops from unquenched Monte Carlo 
simulations at weak couplings, Phys. Rev. D 73 (2006) 094512 [arXiv : hep-lat /0512 012] 



7 



